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ABSTRACT 

We introduce a simple and efficient algorithm for diffusion in smoothed particle hydro- 
dynamics (SPH) simulations and apply it to the problem of chemical mixing. Based 
on the concept of turbulent diffusion, we link the diffusivity of a pollutant to the 
local physical conditions and can thus resolve mixing in space and time. We apply 
our prescription to the evolution of an idealized supernova remnant and find that we 
can model the distribution of heavy elements without having to explicitly resolve hy- 
drodynamic instabilities in the post-shock gas. Instead, the dispersal of the pollutant 
is implicitly modelled through its dependence on the local velocity dispersion. Our 
method can thus be used in any SPH simulation that investigates chemical mixing 
but lacks the necessary resolution on small scales. Potential applications include the 
enrichment of the interstellar medium in present-day galaxies, as well as the inter- 
galactic medium at high redshifts. 



Key words: diffusion 
- early Universe. 



stars: formation - supernova remnants - galaxies: formation 



1 INTRODUCTION 

Understanding chemical enrichment and the dispersal of 
heavy elements in the wake of energetic supernovae (SNe) 
has become essential to a number of fields in astrophysics. 
The details of how enriched material mixes with ambient 
gas are not only relevant for the cooling and fragmenta- 
tion properties of the interstellar medium (ISM), but also 
manifest themselves in the composition and dynamics of 
the resulting stars. As massive stars return metals to the 
ISM, mixing plays a key role in the overall matter cycle 
in galaxies. Tracing this loop back in time, one encoun- 
ters the initial enrichment of the pristine, pure H/He gas 
at the end of the cosmic dark ages, when the very first 
stars, termed Population III (or Pop III), ended their lives 
in violent SN explosions and expelled a significant fraction 



of their mass in metal s dHeger et al. 



120051 : iGreif et all 120071 : IWise fc Abe 



2003; Iwa moto et al.l 



200J). The resulting ume 



2008). The second relates to the enrichment of the inter- 
galac tic medium (IGM) at redshifts z > 5 l|Madau et al.l 
l200ll ). The resulting metallicity distribution is a crucial in- 
gredie nt in modeling the reionization history of the Uni- 
verse l|Yoshida et al.l 12004 iFurlanetto fc Loebl 120051 ) and 
in determining when cosmic star formation shifted from a 
predominantly hi gh-mass, Pop III mode, to a more nor- 
mal Pop II mode jschneider et al.ll2002l:[Mackev et al.ll2003l : 
IVenkatesanll200rj ; rTornatore et al.ll2007l ) . 

Unfortunately, an accurate treatment of mixing, 
while simultaneously simulating the larger-scale environ- 
ment, is presently not feasible, as the turbulent mo- 
tions responsible for mixing typically cascade down to 
very small scales. A frequently encountered approach to 
this problem is to assume that the products of stel- 
lar nucleosynthesis are distributed within a fixed vol- 
e-K- 



i Norman et all 12004 r 



mixing in the early Universe has two aspects. The first 
concerns the enrichment inside the first galaxies, when the 
metals ejected by Pop III stars recollapsed into more mas- 
sive haloes that cooled by atomic hydrogen lines and vig- 
orously mixed with pristin e material in the p r esence of a 
highly turbulent medium (|Wise fc Abell 120071 : IGreif et all 
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.Scannapieco et alj 120051 : 
iKobavashi et al.ll2007l : iTornatore et al.ll2007l ). Significantly 
better results can be achieved in grid-based simulations by 
relating the mass flux between cells to the mixing efficiency, 
even though it remains unclear h ow much of this mi xing is 
numerical instead of physical (e.g. IWise fc Abelllioosl) . Such 
a direct approach is difficult to implement in SPH simula- 
tions, due to their Lagrangian nature, and instead chem- 
ical mixing has been modelle d as a diffusion process (e.g. 
iMartmez-Serrano et al.l 12008) . This is somewhat more ac- 
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curate since the rms displacement of a fluid element in a 
homogeneously and isotropically driven turbulent velocity 
field can be described by the diffusion equation, with the 
diffusion coefficient set by t he velocity dispersion and the 
typical shock travel distance l|Klessen fc Linll2003l ). 

Prescriptions with a constant diffusion coefficient have 
been applie d to models of the chemical evolution of th e 
Milky Way (iKarlssonl 120051 ; iKarlsson fc Gustafssonl 120051) 
galax ies in a cosmological context dMartmez-Serrano et al.l 
|200S|). and the en vironment of the first galaxies 



(Karl sson et al.l 120081 ). In the present paper, we intro- 
duce an improved method that resolves chemical mixing 
in space and time based on the velocity dispersion within 
the SPH smoothing kernel. This yields results in agreement 
with more detailed investigations of mixing in the early 
phases of SN remnants, without having to explicitly resolve 
the hydrodynamic instabilities in the post-shock gas. In 
future work, we plan to use our method to investigate 
the long-term evolution of energetic SNe in a cosmological 
environment, particularly at high redshift when the Uni- 
verse was enriched with the first metals. However, since 
our algorithm is quite generic, we hope that it will prove 
a valuable tool for any SPH simulation that attempts to 
follow the mixing of pollutants. 

The structure of our work is as follows. In Section 2, 
we describe our model for diffusion and its numerical imple- 
mentation in t he cosmological simulation code GADGET-2 
!|Springelll200"5h , followed by a series of idealized test simula- 
tions (Section 3). We then discuss the relation between the 
diffusion coefficient and the velocity dispersion and apply 
our prescription to the evolution of an idealized SN remnant 
(Section 4). Finally, in Section 5 we summarize our results 
and assess their implications. For consistency, all quoted dis- 
tances are physical, unless noted otherwise. 



2 DIFFUSION ALGORITHM 

Diffusion plays an important role in a variety of astro- 
phys ical contexts. Amon g them are thermal conduction 
fe.g.lJubelgas et al.|l20oll) . heat transfer in shear flows (e.g. 



IWadslev et al]|2008h . the microscopic diffusion of particles, 
such as photons in stellar interiors, or the spatial corre- 
lation of individual fluid elements in a turbulent medium 
l|Klessen fc Linll2003l ). These processes are all described by 
the diffusion equation, commonly written in the form: 



dt p y ' 



(1) 



where c is the concentration of a contaminant fluid per unit 
mass, D is the diffusion coefficient, which can be a function 
of space and time, and d/dt the Lagrangian derivative, or the 
derivative following the motion. The diffusion coefficient has 
dimensions M L" 1 T _ , suggesting that it can be represented 
as the product of the local density and some typical length- 
scale and velocity, such as the particle mean free path and 
velocity dispersion in a microscopic picture of diffusion. 



2.1 Numerical implementation 

In the SPH formalism, the diffusion equation can be reduced 
to a discrete summation over all particles within the smooth- 
ing kernel: 

dci 
~dt 



— = > Kijia 



(2) 



with 



Pipj (A + I), 



(3) 



where i and j denote the particle indices, m the mass, p 
the density, Wij the kernel and ry, rij the vector and abso- 
lute separations between part icles i and j, respective ly (for a 
more detailed derivation, see iMonaghan et al.ll2005D . In the 
above equation, the arithmetic mean of the diffusion coef- 
ficient has been replaced by the harmonic mean, which has 
proven to be more robust. Furthermore, the second deriva- 
tive has been replaced by a term involving the gradient and 
the particle separation, since a direct computati on of the 
second derivative is problematic (]Monaghanll2005h . 

The solution to equation (2) can either be determined 
explicitly, which requires an additional constraint on the 
time-step to ensure numerical stability, or implicitly, which 
requires the solution of a coupled set of differential equations 
involving all 'active' particles (i.e. all SPH particles that are 
being updated on the current time-step). The explicit ap- 
proach is extremely difficult to implement in a conservative 
fashion in an SPH code that allows individual particles to 
have different time-steps (e.g. GADGET-2). This is because 
in a code of this type, neighbouring particles will sometimes 
have different time-steps. Consequently, the increase in con- 
centration at particle i caused by diffusion from particle j 
will sometimes be computed at a different time from the 
corresponding decrease in concentration at particle j caused 
by diffusion to particle i. To ensure conservation, the in- 
crease and decrease must exactly balance, but in general 
they will not if i and j have different time-steps. One could, 
of course, avoid this problem by ensuring that all particles 
are synchronized before their concentrations are updated. 
However, the required synchronizations would have to oc- 
cur very frequently, and so one would lose essentially all of 
the benefits gained by allowing the particles to have individ- 
ual time-steps. A further undesirable feature of the explicit 
approach is the fact that the time-steps required to stably 
model the diffusion can become very small. Consideration of 
equation (I) shows that the required time-step scales with 
the spatial resolution - represented in an SPH code by the 
smoothing length h - as 



At oc h 



(4) 



In comparison, the standard Courant time-step scales only 
linearly with h. 

An implicit approach to the solution of equation (2) 
avoids some of these problems, as it allows one to take larger 
time-steps without compromising numerical stability. How- 
ever, this comes at a cost: the coupled differential equations 
representing the diffusion must be solved iteratively, and it 
is difficult to do this in a fashion that can be efficiently 
parallelized. In addition, one still has to deal with the syn- 
chronization problem discussed above. 



© 0000 RAS, MNRAS 000, 000-000 



Chemical mixing in SPH simulations 3 



In view of the disadvantages of both standard ap- 
proaches, it is interesting to explore the viability of simpler, 
but more approximate approaches, such as the one presented 
in this paper. To obtain our approximation, we assume that 
the densities and concentrations of all active particles do 
not change significantly over a time interval At, allowing a 
direct integration of equation (2): 

B . 



Ci(t + At) = Ci{t )e 



with 



AAt 



)> 



j 

and 

B = J2K i: 



(5) 



(6) 



(7) 



For large time-steps, the concentration of particle i thus 
tends to the average among its neighbours, while for small 
time-steps it remains close to its original concentration. 

We implement this approach by performing the required 
operations at a global synchronization point in the density 
routine of GADGET-2. After the new densities and smooth- 
ing lengths have been computed, we update the coefficients 
Kij and subsequently use equation (5) to determine the new 
concentrations of all active SPH particles. In a final step, 
we renormalize the obtained concentrations with a global 
factor such that the total concentration is conserved. Since 
the Courant condition does not allow for significant changes 
in density between time-steps, and diffusion generally pro- 
gresses slower than the speed of sound, our implementation 
is quite generic and can be applied to a number of prob- 
lems in astrophysics. We have found that the algorithm is 
remarkably stable even for very short diffusion timescales, 
since particles tend to equilibrate their concentrations and 
neighbouring particles are generally active at the same time. 
The additional CPU consumption is minimal since we uti- 
lize the pre-existing neighbour search. By the same token, 
the algorithm is easy to implement in any SPH code and is 
not restricted to GADGET-2. 



2.2 Test problems 

In this section we investigate the formal accuracy of the algo- 
rithm by performing a number of idealized test simulations. 
We initialize all simulations in a periodic, uniform density 
box with 1 million SPH particles and length and sound- 
crossing time set to unity, such that we may conveniently 
quote the elapsed time in units of the sound-crossing time. 
We adopt a hydrogen number density of «,h = 1 cm -3 and 
a mean molecular weight corresponding to that of a neutral, 
primordial gas. We place the particles on a grid with a very 
small random displacement, and in the first test problem 
also consider a fully random distribution of particles, such 
that the density fluctuates considerably around the mean. 
This gives a better handle on the performance of our algo- 
rithm under more realistic circumstances. In all cases, we 
use the maximum diffusivity that is accurately modelled by 
our algorithm, determined by the Courant condition (see 
Sections 2.1 and 3.1): 



where c s is the sound speed and p and h are determined 
by the mean density tih = 1 cm~ 3 , such that the diffusion 
coefficient becomes a fixed numerical value. 

In the first test problem, we consider the propagation 
of an initial (^-distribution, to which the analytic solution is 
Green's function of the diffusion equation: 



G(x,x',t) 



(27TCT 2 ) 3 / 2 



exp ( ) (9) 



with variance 



(10) 

This configuration is reproduced by setting the concentra- 
tion of the central particle to unity, and of all others to zero. 
In the left-hand column of Fig. 1, we compare the analytic 
solution to the simulation results at three different output 
times. Early on, diffusion progresses somewhat too rapidly, 
since the concentrations of neighbouring particles differ sub- 
stantially. In this case, equation (5) effectively breaks down, 
but the resulting deviations remain small and do not influ- 
ence the long-term behavior, where all three curves become 
nearly indistinguishable. Note that the random distribution 
performs almost as well as the grid-based distribution, show- 
ing that the diffusivity remains unchanged even for high 
density fluctuations. 

The second test problem consists of two individual 5- 
distributions initially separated by Ax = 1/3. The analytic 
solution can be obtained by a convolution of Green's func- 
tion with the initial state of the system: 



c(x,t + At) = J G(x,x',t + At) c(x ,to) dx , 



(11) 



such that the general solution is given by the superposition 
of two Gaussian distributions centred at x = 1/3 and 2/3. 
Similarly, the solution to the third test problem, consisting 
of a slab of uniform concentration between x = 1/3 and 2/3, 
is given by: 



c(x, t) = i erfc ( 
for x < 1/2 and 
c(x, t) = i erfc ( - ■ 



(12) 



(13) 



D = p ft c s 



(8) 



for x > 1/2. In the middle and right-hand columns of Fig. 1, 
we compare the simulation results to the analytic solution, 
showing that there are only minor deviations at early times, 
similar in nature to those found in the first test problem. A 
slice through all three boxes (Fig. 2) shows the solution in 
two dimensions. 

Finally, in a series of resolution studies performed with 
32 3 , 64 3 and 128 3 particles, we have found no correlation 
between the resolution and the deviation from the analytic 
solution. In fact, in all cases, the deviations were compara- 
ble to those found in previous test problems. Considering 
its formal simplicity, the algorithm thus performs remark- 
ably well. Although the diffusivity is initially slightly over- 
predicted for cases in which the diffusion time is shorter 
than the sound-crossing time, we nevertheless find a correct 
long-term behavior. 
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Figure 1. The propagation of a single 5-distribution (left-hand column), two ^-distributions (middle column) and a slab of uniform 
concentration (right-hand column) into an otherwise pristine medium, shown for a grid-based particle distribution (solid lines) and a 
random particle distribution (dashed line), compared to the analytic solution (dotted lines). The temporal evolution is depicted from 
top to bottom and quoted in units of the sound-crossing time. In all three cases, the simulations reproduce the analytic solution, except 
at early times when the underlying assumption of constant concentration for neighbouring particles is not well justified (sec Section 2). 
However, as is evident from the bottom two rows, this does not affect the long-term behavior. 



3 APPLICATION TO CHEMICAL MIXING 

In the previous section, we introduced an SPH formalism 
for diffusion that can be applied to most problems that are 
governed by a diffusion equation. In this section, we focus on 
an application that is particularly relevant to astrophysics: 
the mixing of chemical elements. 

3.1 Chemical mixing as turbulent diffusion 

As a first step towards a model for chemical mixing, one 
must find a connection between t he diffusivity of a poll utant 
and the local physical conditions. iKlessen fc Lml (|2003l ) have 



provided this link by showing that the probability of finding 
a parcel of gas at a given location in a homogeneously and 
isotropically driven turbulent velocity field can be described 
by the diffusion equation, with the diffusion coefficient set 
by the velocity dispersion v and the turbulent driving length 
I: 

D = 2pvl. (14) 

This corresponds to the classical mixing-length approach ex- 
tended into the supersonic regime. If one replaces the proba- 
bility distribution with the concentration of a pollutant, this 
formalism can be reinterpreted to describe chemical mixing. 
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Test Case (b) 



Test Case (c) 
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Figure 2. The propagation of a single <5-distribution (left-hand column), two 5-distributions (middle column) and a slab of uniform 
concentration (right-hand column) into an otherwise pristine medium, shown at a representative output time. 



The sole remaining task is then to provide the parameters 
v and I to the diffusion algorithm, such that the diffusion 
coefficient can become a function of space and time. 

To obtain the necessary parameters, we determine the 
rms velocity dispersion for particle i within its smoothing 
length: 



1 



N, 



ngb 



E 

3 



V, 



(15) 



where iV n gb is the number of neighbours and V; and Vj are 
the bulk velocities of particles i and j. Finally, we equate the 
turbulent driving length with the smoothing length, since 
this is the minimum scale where turbulent motions can be 
resolved. This then yields for the diffusion coefficient: 



Di — 2 pi Vi hi 



(16) 



As desired, the efficiency of chemical mixing is now governed 
entirely by local physical quantities. 

An important underlying assumption of this method is 
that the velocity field on subresolution scales corresponds to 
a homogeneously and isotropically driven turbulent medium, 
i.e. only the magnitude of the velocity dispersion and not its 
three-dimensional structure is taken into account. A further 
implicit assumption is that the mixing efficiency on sub- 
resolution scales is set by the corresponding value on the 
scale of the smoothing kernel. For these reasons our method 
yields an upper limit to the amount of mixing that can oc- 
cur, and does not capture the details of a real turbulent 
cascade. However, since turbulent motions on scales larger 
than the smoothing length are explicitly resolved, our ap- 
proach should suffice for most practical purposes. 

We implement the above steps in our algorithm by per- 
forming a previous neighbour search that finds the velocity 
dispersion for all active SPH particles, which is then used in 
equation (3) to obtain the coefficients Kij. This closes the 
required set of equations, and in the following subsection we 
use our complete model to investigate the mixing of metals 
in a SN remnant. 



3.2 Mixing in supernova remnants 

The mixing of gas in the post-shock regions of SN remnants 
has previously been investigated, leading to the consen- 
sus view that secondary shocks trigger Rayleigh- Taylor and 
Kelvin-Helmholtz instabilities that mix pristine g as from the 
dense shell with alre ady enriched material (e.g. lGulllll973l ; 
IChevalier et al.lll992n . The aim of our algorithm is to cap- 
ture this mixing without having to explicity resolve the cor- 
responding hydrodynamic i nstabilities, which cannot easil y 
be modelled using SPH (e.g. lAgertz et alj|2007l ; |Pricell2008l ) . 
To verify its ability to do this, we perform a test simulation 
of an idealized SN explosion with the same setup as in Sec- 
tion 2.2. We distribute an explosion energy of 10 51 erg as 
thermal energy to the iV ng b nearest neighbours around the 
centre of the box, which results in the formation of a shock 
and reproduces the ideal Sedo v- Taylor solution of a blast 
wave in a uniform medium (see lGreif et al. 2007). Further- 
more, we set the initial metallicities of the central particles 
to unity. 

In Fig. 3, we show the radial velocity, diffusion coeffi- 
cient and metallicity as a function of distance from the SN 
progenitor at three different output times. The high veloc- 
ity dispersion near the forward shock is responsible for the 
elevated diffusion coefficient, which leads to efficient chemi- 
cal mixing of pristine gas passing through the shock. As is 
evident from Fig. 3, most parcels of gas behind the shock 
have nearly equilibrated their metallicities. Our algorithm 
thus implicitly accounts for mixing due to the dependence 
of the diffusion coefficient on the local velocity dispersion. 
Intriguingly, we also find metals ahead of the forward shock. 
Although this initially seems unphysical, it should be re- 
membered that our algorithm aims to capture mixing caused 
by unresolved instabilities, and that in a real SN remnant, 
it is these instabilities that would transport metals ahead of 
the mean position of the shock. Even further out, the degree 
of enrichment drops roughly exponentially, since the velocity 
dispersion, and hence the diffusion coefficient, tend to zero. 
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Velocity Diffusion coefficient Metallicity 




Figure 3. The radial velocity (left-hand column), diffusion coefficient (middle column) and metallicity (right-hand column) as a function 
of distance from the SN progenitor, with the position of the forward shock denoted by the dotted line. The temporal evolution is depicted 
from top to bottom and quoted in units of the sound-crossing time. The high velocity dispersion near the main shock leads to an elevated 
diffusion coefficient and thus very efficient chemical mixing. Due to the dependence on velocity dispersion, mixing by hydrodynamic 
instabilities in the remnant is implicitly accounted for. 



A second mechanism becomes important once the SN 
remnant has stalled: gas from the dense shell ahead of the 
shock falls back on to the remnant, becomes Rayleigh- Taylor 
unstable, and vigorously mixes with the interior gas. This is 
especially effective in a realistic cosmological setting, where 
the ce ntral potential we l l deepens as the S N remnant ex- 
pands (|Greif et all 120071 : IWise fc Abellbooih . However, we 
postpone a detailed treatment of this issue to dedicated 
high-resolution simulations in future work. 



4 SUMMARY AND CONCLUSIONS 

We have introduced a simple and efficient algorithm for dif- 
fusion in SPH simulations and investigated its accuracy with 
a number of idealized test cases. Although the algorithm is 
quite generic and can be applied to most problems involving 
diffusion, we have here focused on a model for the dispersal 
of enriched material in supernova ex plosions. Adopt i ng th e 
mixing length approach discussed bv lKlessen fc LirJ |2003), 
we link the diffusivity of the pollutant to the local phys- 
ical conditions and can thus describe the space- and time- 
dependent mixing process. We have applied our prescription 
to an idealized SN explosion and found that we can properly 
describe the mixing process without having to explicitly re- 
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solve the hydro-dynamic instabilities in the post-shock gas. 
Instead, this process is implicitly governed by the depen- 
dence on the local velocity dispersion. Our method can thus 
be used in any SPH simulation that attempts to follow the 
mixing of a pollutant but lacks the necessary resolution on 
small scales. This will be relevant, in particular, for simula- 
tions that aim at simultaneously representing the large-scale 
environment of the SN explosion and the fine-grained mix- 
ing. A crucial assumption of our model is that resolved mo- 
tions cascade down to subresolution scales which then homo- 
geneously mix the gas. For this reason our method yields an 
upper limit to the mixing efficiency and the resulting degree 
of homogeneity, and cannot capture the details of a real tur- 
bulent cascade. However, since turbulent motions on scales 
larger than the smoothing length are explicitly resolved, this 
approximation should generally be valid. 

Our new algorithm ties in well with the current gen- 
eration of radiation-hydrodynamical simulations of star for- 
mation, both at high redshifts and in the local Universe. 
The ultimate goal is to describe the interrelated cycle of gas 
fragmentation and stellar feedback. To achieve adequate re- 
alism, radiative effects have to be considered together with 
the mechanical and chemical feedback due to SN explosions. 
The methodology developed here allows us to include chem- 
ical feedback, suitable even for extremely large simulations. 
We will report on specific applications in future studies. 
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